options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

non linear regression with one explanatory variable

fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

ex8-0.stan

single linear regression

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)
## Initial log joint probability = -275175 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       36      -38.1623    0.00367238   0.000290015           1           1       86    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -38.16
##    b0         7.96
##    b1         1.80
##    s          4.09
##    m0[1]     33.49
##    m0[2]     37.42
##    m0[3]     25.27
##    m0[4]     43.10
##    m0[5]     36.89
##    m0[6]     25.44
##    m0[7]     41.35
##    m0[8]     16.79
##    m0[9]     29.90
##    m0[10]    19.98
##    m0[11]    25.63
##    m0[12]    29.00
##    m0[13]    38.77
##    m0[14]    37.31
##    m0[15]    39.65
##    m0[16]    32.01
##    m0[17]    34.18
##    m0[18]    13.78
##    m0[19]    22.34
##    m0[20]    22.12
##    y0[1]     30.23
##    y0[2]     37.79
##    y0[3]     32.51
##    y0[4]     41.73
##    y0[5]     44.44
##    y0[6]     20.67
##    y0[7]     38.04
##    y0[8]     16.54
##    y0[9]     31.74
##    y0[10]    19.97
##    y0[11]    24.49
##    y0[12]    35.25
##    y0[13]    28.08
##    y0[14]    27.72
##    y0[15]    44.25
##    y0[16]    31.67
##    y0[17]    38.49
##    y0[18]    14.84
##    y0[19]    10.39
##    y0[20]    22.31
##    m1[1]     13.78
##    m1[2]     17.04
##    m1[3]     20.29
##    m1[4]     23.55
##    m1[5]     26.81
##    m1[6]     30.07
##    m1[7]     33.33
##    m1[8]     36.58
##    m1[9]     39.84
##    m1[10]    43.10
##    y1[1]      9.63
##    y1[2]     16.29
##    y1[3]     22.19
##    y1[4]     22.72
##    y1[5]     27.07
##    y1[6]     30.94
##    y1[7]     26.22
##    y1[8]     31.60
##    y1[9]     42.91
##    y1[10]    42.25
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -38.26 -37.96 1.22 1.04 -40.65 -36.95 1.01      502      557
##    b0       8.36   8.22 2.87 2.83   3.74  13.12 1.01      317      316
##    b1       1.78   1.78 0.22 0.21   1.40   2.13 1.01      348      396
##    s        4.59   4.50 0.81 0.74   3.50   6.11 1.00      945      911
##    m0[1]   33.50  33.51 1.10 1.04  31.71  35.30 1.01     2024     1283
##    m0[2]   37.37  37.37 1.34 1.28  35.17  39.51 1.01     1403     1221
##    m0[3]   25.41  25.43 1.18 1.14  23.51  27.33 1.01      510      680
##    m0[4]   42.97  42.96 1.86 1.75  39.86  45.97 1.00      713      975
##    m0[5]   36.85  36.87 1.30 1.24  34.72  38.94 1.01     1509     1205
##    m0[6]   25.58  25.59 1.17 1.14  23.69  27.48 1.00      522      695
##    m0[7]   41.25  41.23 1.69 1.61  38.46  43.93 1.00      818      976
##    m0[8]   17.05  16.99 1.91 1.89  13.94  20.24 1.01      333      357
##    m0[9]   29.97  29.97 1.03 0.99  28.29  31.68 1.01     1300     1196
##    m0[10]  20.20  20.17 1.60 1.58  17.62  22.83 1.01      357      398
##    m0[11]  25.76  25.78 1.16 1.13  23.90  27.64 1.00      536      704
##    m0[12]  29.08  29.08 1.04 1.00  27.40  30.78 1.01     1055     1156
##    m0[13]  38.71  38.71 1.46 1.38  36.31  41.02 1.00     1180     1120
##    m0[14]  37.26  37.26 1.34 1.27  35.08  39.39 1.01     1424     1172
##    m0[15]  39.57  39.57 1.53 1.46  37.04  42.01 1.00     1009     1121
##    m0[16]  32.04  32.04 1.05 0.99  30.35  33.75 1.01     1872     1267
##    m0[17]  34.19  34.18 1.13 1.06  32.33  36.05 1.01     1988     1286
##    m0[18]  14.09  13.98 2.23 2.19  10.46  17.74 1.01      323      364
##    m0[19]  22.52  22.50 1.40 1.38  20.28  24.79 1.01      396      445
##    m0[20]  22.30  22.27 1.41 1.40  20.03  24.62 1.01      391      445
##    y0[1]   33.65  33.62 4.83 4.58  25.71  41.79 1.00     1877     1727
##    y0[2]   37.27  37.32 4.82 4.64  29.31  45.20 1.00     2025     1864
##    y0[3]   25.39  25.46 4.76 4.61  17.28  32.85 1.00     1616     1816
##    y0[4]   43.04  43.15 5.00 4.71  34.87  51.06 1.00     1836     1970
##    y0[5]   37.06  37.09 4.77 4.74  29.41  44.75 1.00     1965     1835
##    y0[6]   25.75  25.81 4.77 4.71  17.94  33.54 1.00     1833     1947
##    y0[7]   41.21  41.16 5.05 4.93  32.96  49.55 1.00     1689     1982
##    y0[8]   16.94  17.02 5.00 4.86   8.80  24.96 1.00     1343     1548
##    y0[9]   30.06  30.16 4.67 4.67  22.39  37.84 1.00     2204     1976
##    y0[10]  20.23  20.35 4.91 4.60  11.60  28.22 1.00     1341     1648
##    y0[11]  25.66  25.65 4.79 4.82  17.88  33.49 1.00     1906     1866
##    y0[12]  29.02  29.09 4.60 4.65  21.33  36.48 1.00     1963     1743
##    y0[13]  38.79  38.54 4.88 4.50  31.21  46.90 1.00     1788     2001
##    y0[14]  37.17  37.27 4.72 4.54  29.27  44.77 1.00     1983     2124
##    y0[15]  39.51  39.61 4.87 4.73  31.25  47.33 1.00     1794     1876
##    y0[16]  32.17  32.18 4.84 4.73  24.18  39.96 1.00     1951     1921
##    y0[17]  34.23  34.17 4.90 4.77  26.04  42.18 1.00     2098     1878
##    y0[18]  14.03  14.01 5.15 4.79   5.57  22.73 1.00     1001     1218
##    y0[19]  22.68  22.81 4.86 4.57  14.75  30.60 1.00     1377     1521
##    y0[20]  22.31  22.27 4.77 4.60  14.82  30.09 1.00     1620     1743
##    m1[1]   14.09  13.98 2.23 2.19  10.46  17.74 1.01      323      364
##    m1[2]   17.29  17.24 1.89 1.87  14.23  20.43 1.01      334      371
##    m1[3]   20.50  20.47 1.57 1.55  17.98  23.10 1.01      361      417
##    m1[4]   23.71  23.70 1.30 1.25  21.62  25.83 1.01      430      485
##    m1[5]   26.92  26.92 1.10 1.10  25.14  28.71 1.00      644      879
##    m1[6]   30.13  30.13 1.03 0.99  28.46  31.84 1.01     1352     1174
##    m1[7]   33.34  33.34 1.09 1.03  31.56  35.12 1.01     2020     1265
##    m1[8]   36.55  36.56 1.28 1.22  34.42  38.61 1.01     1579     1259
##    m1[9]   39.76  39.75 1.55 1.48  37.23  42.22 1.00      979     1121
##    m1[10]  42.97  42.96 1.86 1.75  39.86  45.97 1.00      713      975
##    y1[1]   14.18  14.13 5.13 4.95   5.80  22.44 1.00     1036     1540
##    y1[2]   17.32  17.44 5.12 4.83   8.75  25.43 1.00     1160     1717
##    y1[3]   20.57  20.63 4.83 4.67  12.64  28.35 1.00     1703     1752
##    y1[4]   23.70  23.81 4.79 4.60  15.65  31.42 1.00     1669     1799
##    y1[5]   26.75  26.72 4.84 4.69  18.99  34.83 1.00     1998     1808
##    y1[6]   30.12  30.05 4.77 4.70  22.49  38.06 1.00     1811     1931
##    y1[7]   33.47  33.53 4.71 4.64  25.76  40.97 1.00     2046     1808
##    y1[8]   36.50  36.49 4.95 4.62  28.28  44.60 1.00     2008     1922
##    y1[9]   39.80  39.82 4.87 4.52  31.71  47.78 1.00     1886     1972
##    y1[10]  42.84  42.74 5.02 4.93  34.84  51.31 1.00     1820     1738

quadratic regression

y=b0+b2(x-b1)**2

ex8-1.stan

quadratic regression

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real b2;
  real<lower=0> s;
}
model{
  y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b2*(x[i]-b1)^2;
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b2*(x1[i]-b1)^2;
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)
## Initial log joint probability = -574452 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpgKfISJ/model-14463151a723.stan', line 15, column 2 to column 29) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpgKfISJ/model-14463151a723.stan', line 15, column 2 to column 29) 
## Error evaluating model log probability: Non-finite gradient. 
##       72      -8.59084   0.000254104    0.00246273           1           1      152    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -8.59
##    b0         5.11
##    b1         4.01
##    b2         0.52
##    s          0.93
##    m0[1]      5.64
##    m0[2]     10.65
##    m0[3]     17.22
##    m0[4]      5.14
##    m0[5]      7.61
##    m0[6]      8.58
##    m0[7]     10.02
##    m0[8]     14.83
##    m0[9]      5.51
##    m0[10]     5.12
##    m0[11]     5.46
##    m0[12]     5.20
##    m0[13]    17.65
##    m0[14]     9.01
##    m0[15]     6.50
##    m0[16]    12.46
##    m0[17]     5.94
##    m0[18]     7.64
##    m0[19]     7.53
##    m0[20]     7.08
##    y0[1]      5.18
##    y0[2]     10.08
##    y0[3]     18.66
##    y0[4]      6.64
##    y0[5]      8.46
##    y0[6]      8.90
##    y0[7]     10.15
##    y0[8]     15.35
##    y0[9]      5.81
##    y0[10]     4.79
##    y0[11]     4.38
##    y0[12]     3.22
##    y0[13]    18.89
##    y0[14]     8.11
##    y0[15]     6.22
##    y0[16]    13.23
##    y0[17]     5.67
##    y0[18]     5.87
##    y0[19]     9.43
##    y0[20]     7.78
##    m1[1]      8.58
##    m1[2]      6.70
##    m1[3]      5.55
##    m1[4]      5.11
##    m1[5]      5.40
##    m1[6]      6.40
##    m1[7]      8.13
##    m1[8]     10.58
##    m1[9]     13.75
##    m1[10]    17.65
##    y1[1]      6.38
##    y1[2]      6.04
##    y1[3]      5.65
##    y1[4]      6.04
##    y1[5]      6.12
##    y1[6]      6.18
##    y1[7]      8.22
##    y1[8]     10.12
##    y1[9]     13.70
##    y1[10]    15.84
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -10.91 -10.52 1.56 1.38 -13.94 -9.07 1.00      616      918
##    b0       5.11   5.12 0.37 0.36   4.47  5.70 1.01      459      553
##    b1       3.99   4.00 0.15 0.14   3.72  4.22 1.00     1056      870
##    b2       0.52   0.51 0.05 0.05   0.43  0.60 1.00      750     1058
##    s        1.10   1.06 0.20 0.18   0.82  1.47 1.00      763      950
##    m0[1]    5.64   5.64 0.33 0.32   5.10  6.18 1.01      537      748
##    m0[2]   10.66  10.65 0.37 0.35  10.07 11.27 1.00     1201     1220
##    m0[3]   17.19  17.20 0.63 0.58  16.13 18.20 1.00     1635     1267
##    m0[4]    5.16   5.17 0.38 0.37   4.49  5.77 1.01      473      514
##    m0[5]    7.58   7.58 0.47 0.45   6.81  8.36 1.00     1186      754
##    m0[6]    8.54   8.54 0.58 0.57   7.59  9.48 1.00     1259      694
##    m0[7]   10.03  10.03 0.36 0.35   9.44 10.63 1.00     1015     1341
##    m0[8]   14.81  14.82 0.50 0.46  13.95 15.64 1.00     2130     1319
##    m0[9]    5.53   5.55 0.40 0.39   4.83  6.16 1.01      502      525
##    m0[10]   5.14   5.14 0.36 0.35   4.50  5.71 1.01      464      606
##    m0[11]   5.48   5.50 0.40 0.38   4.79  6.11 1.01      499      525
##    m0[12]   5.23   5.24 0.39 0.37   4.54  5.85 1.01      478      529
##    m0[13]  17.61  17.62 0.66 0.61  16.51 18.68 1.00     1543     1233
##    m0[14]   9.03   9.03 0.37 0.36   8.42  9.63 1.00      815      918
##    m0[15]   6.49   6.49 0.37 0.36   5.87  7.08 1.00      854      686
##    m0[16]  12.46  12.45 0.40 0.37  11.78 13.13 1.01     1918     1337
##    m0[17]   5.94   5.94 0.34 0.33   5.38  6.48 1.01      619      771
##    m0[18]   7.66   7.67 0.38 0.38   7.03  8.27 1.00      649      780
##    m0[19]   7.50   7.50 0.46 0.44   6.74  8.27 1.00     1174      720
##    m0[20]   7.06   7.07 0.42 0.40   6.36  7.75 1.00     1078      694
##    y0[1]    5.62   5.66 1.16 1.09   3.74  7.52 1.00     1635     1730
##    y0[2]   10.62  10.66 1.15 1.13   8.64 12.42 1.00     1972     1809
##    y0[3]   17.21  17.22 1.33 1.33  15.05 19.33 1.00     1940     1823
##    y0[4]    5.20   5.21 1.20 1.12   3.20  7.10 1.00     1597     1449
##    y0[5]    7.60   7.65 1.23 1.13   5.52  9.57 1.00     1785     1535
##    y0[6]    8.52   8.51 1.24 1.25   6.47 10.57 1.00     1682     1698
##    y0[7]   10.04  10.08 1.21 1.15   8.00 12.05 1.00     1851     1772
##    y0[8]   14.83  14.81 1.22 1.15  12.79 16.85 1.00     1957     1747
##    y0[9]    5.54   5.53 1.17 1.12   3.60  7.41 1.00     1591     1743
##    y0[10]   5.14   5.14 1.12 1.05   3.25  6.92 1.00     1473     1775
##    y0[11]   5.49   5.52 1.14 1.10   3.61  7.34 1.00     1618     1531
##    y0[12]   5.22   5.23 1.16 1.09   3.33  7.12 1.00     1656     1805
##    y0[13]  17.62  17.60 1.30 1.22  15.47 19.72 1.00     2038     1909
##    y0[14]   9.05   9.06 1.20 1.14   7.10 11.00 1.00     1788     1788
##    y0[15]   6.50   6.52 1.15 1.10   4.58  8.37 1.00     2026     1845
##    y0[16]  12.45  12.47 1.19 1.12  10.54 14.46 1.00     1993     1601
##    y0[17]   5.94   5.92 1.13 1.11   4.06  7.77 1.00     1822     1909
##    y0[18]   7.68   7.70 1.17 1.10   5.78  9.59 1.00     1801     1721
##    y0[19]   7.54   7.51 1.18 1.16   5.59  9.45 1.00     1641     1912
##    y0[20]   7.09   7.09 1.18 1.14   5.20  9.06 1.00     1956     1821
##    m1[1]    8.54   8.54 0.58 0.57   7.59  9.48 1.00     1259      694
##    m1[2]    6.69   6.69 0.38 0.37   6.04  7.31 1.00      948      682
##    m1[3]    5.55   5.55 0.33 0.32   5.01  6.08 1.01      517      734
##    m1[4]    5.13   5.13 0.36 0.35   4.49  5.71 1.01      465      589
##    m1[5]    5.42   5.44 0.40 0.38   4.73  6.05 1.01      494      531
##    m1[6]    6.43   6.45 0.40 0.39   5.74  7.05 1.01      557      604
##    m1[7]    8.15   8.16 0.38 0.37   7.52  8.76 1.00      699      675
##    m1[8]   10.59  10.59 0.37 0.35  10.00 11.20 1.00     1177     1218
##    m1[9]   13.74  13.75 0.45 0.41  12.96 14.47 1.00     2164     1303
##    m1[10]  17.61  17.62 0.66 0.61  16.51 18.68 1.00     1543     1233
##    y1[1]    8.54   8.54 1.28 1.27   6.45 10.62 1.00     1706     1761
##    y1[2]    6.67   6.70 1.19 1.11   4.64  8.54 1.00     1827     1962
##    y1[3]    5.56   5.54 1.15 1.09   3.67  7.42 1.00     1598     1520
##    y1[4]    5.13   5.16 1.18 1.11   3.25  7.07 1.00     1696     1658
##    y1[5]    5.39   5.39 1.16 1.09   3.51  7.31 1.00     1695     1722
##    y1[6]    6.38   6.37 1.17 1.14   4.42  8.24 1.00     1746     1580
##    y1[7]    8.15   8.16 1.18 1.16   6.16  9.96 1.00     1642     1490
##    y1[8]   10.61  10.61 1.17 1.13   8.71 12.53 1.00     1965     1902
##    y1[9]   13.71  13.69 1.23 1.18  11.67 15.70 1.00     1986     1863
##    y1[10]  17.58  17.58 1.31 1.21  15.51 19.76 1.00     1938     1841

both log regression

log y=b0+b1log x x,y>0 y=exp b0 x**b1

ex8-2.stan

both log regression

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector<lower=0>[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~lognormal(b0+b1*log(x),s);
}
generated quantities{
  vector[N] m0;
  vector<lower=0>[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*log(x[i]);
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector<lower=0>[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*log(x1[i]);
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
             qplot(log(x),log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -176.919 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22      -24.7469   0.000170539   2.69434e-06           1           1       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -24.75
##    b0         1.00
##    b1         1.84
##    s          0.83
##    m0[1]      2.17
##    m0[2]      4.86
##    m0[3]      3.00
##    m0[4]      2.61
##    m0[5]      3.65
##    m0[6]      5.07
##    m0[7]      4.83
##    m0[8]      2.22
##    m0[9]     -0.85
##    m0[10]     1.56
##    m0[11]     1.94
##    m0[12]     5.21
##    m0[13]     0.94
##    m0[14]     0.16
##    m0[15]     4.94
##    m0[16]    -0.11
##    m0[17]    -2.66
##    m0[18]     4.51
##    m0[19]     4.26
##    m0[20]     4.57
##    y0[1]      5.03
##    y0[2]    114.89
##    y0[3]     16.36
##    y0[4]      9.55
##    y0[5]     42.82
##    y0[6]    301.65
##    y0[7]    366.15
##    y0[8]      9.87
##    y0[9]      0.75
##    y0[10]     5.68
##    y0[11]     2.39
##    y0[12]   173.30
##    y0[13]     1.52
##    y0[14]     0.25
##    y0[15]    85.51
##    y0[16]     0.53
##    y0[17]     0.09
##    y0[18]    42.11
##    y0[19]   516.22
##    y0[20]    89.65
##    m1[1]     -2.66
##    m1[2]      1.36
##    m1[3]      2.53
##    m1[4]      3.24
##    m1[5]      3.75
##    m1[6]      4.15
##    m1[7]      4.48
##    m1[8]      4.76
##    m1[9]      5.00
##    m1[10]     5.21
##    y1[1]      0.08
##    y1[2]      2.82
##    y1[3]     12.47
##    y1[4]     61.17
##    y1[5]     70.81
##    y1[6]     24.59
##    y1[7]     54.41
##    y1[8]    335.42
##    y1[9]    232.40
##    y1[10]   382.55
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median     sd    mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -26.47 -26.12   1.30   1.05 -28.90 -25.10 1.00      904     1011
##    b0       1.01   1.01   0.26   0.25   0.59   1.41 1.01      823     1046
##    b1       1.83   1.83   0.17   0.17   1.55   2.11 1.00     1133     1191
##    s        0.94   0.91   0.17   0.15   0.70   1.27 1.00     1260     1200
##    m0[1]    2.18   2.17   0.21   0.21   1.82   2.51 1.00     1030     1081
##    m0[2]    4.86   4.85   0.30   0.28   4.35   5.36 1.00     2293     1652
##    m0[3]    3.00   3.00   0.21   0.20   2.65   3.35 1.00     1473     1465
##    m0[4]    2.62   2.62   0.21   0.19   2.27   2.95 1.00     1243     1297
##    m0[5]    3.65   3.65   0.23   0.21   3.25   4.04 1.00     1970     1438
##    m0[6]    5.06   5.06   0.31   0.30   4.54   5.59 1.00     2264     1596
##    m0[7]    4.83   4.83   0.29   0.28   4.33   5.33 1.00     2298     1652
##    m0[8]    2.22   2.22   0.21   0.21   1.87   2.56 1.00     1059     1078
##    m0[9]   -0.84  -0.84   0.39   0.37  -1.46  -0.21 1.01      821      889
##    m0[10]   1.56   1.56   0.23   0.23   1.18   1.92 1.01      878      997
##    m0[11]   1.94   1.94   0.22   0.21   1.58   2.29 1.01      957     1110
##    m0[12]   5.21   5.20   0.32   0.31   4.66   5.75 1.00     2233     1546
##    m0[13]   0.95   0.95   0.26   0.25   0.52   1.36 1.01      820     1002
##    m0[14]   0.17   0.17   0.31   0.29  -0.34   0.68 1.01      803      869
##    m0[15]   4.93   4.93   0.30   0.29   4.42   5.45 1.00     2285     1681
##    m0[16]  -0.09  -0.09   0.33   0.31  -0.63   0.45 1.01      805      869
##    m0[17]  -2.64  -2.64   0.53   0.53  -3.50  -1.75 1.01      872      888
##    m0[18]   4.50   4.50   0.27   0.26   4.04   4.97 1.00     2325     1596
##    m0[19]   4.26   4.25   0.26   0.24   3.82   4.70 1.00     2296     1651
##    m0[20]   4.57   4.57   0.28   0.26   4.10   5.05 1.00     2320     1596
##    y0[1]   14.53   8.80  19.98   7.58   1.87  43.99 1.00     1984     1892
##    y0[2]  205.52 129.33 254.86 112.09  25.23 630.10 1.00     1802     1800
##    y0[3]   32.26  20.22  49.33  16.91   4.18  89.64 1.00     2058     1859
##    y0[4]   22.71  13.34  32.73  10.89   2.82  68.87 1.00     1871     1892
##    y0[5]   60.24  37.58  78.61  31.25   7.61 173.87 1.00     2165     2100
##    y0[6]  273.64 153.51 434.34 135.14  31.94 846.11 1.00     2005     1800
##    y0[7]  210.84 123.07 368.38 104.13  25.80 627.61 1.00     1957     1877
##    y0[8]   14.64   9.24  20.65   7.90   1.78  44.43 1.00     1948     1778
##    y0[9]    0.76   0.42   1.76   0.38   0.08   2.33 1.00     1593     1945
##    y0[10]   7.83   4.80  13.55   4.05   0.96  22.34 1.00     1562     1965
##    y0[11]  11.47   7.04  18.32   5.94   1.45  32.76 1.00     2032     1769
##    y0[12] 327.15 188.81 685.13 163.01  33.92 968.30 1.00     2096     1973
##    y0[13]   4.53   2.73   6.93   2.37   0.47  13.75 1.00     1868     1723
##    y0[14]   2.06   1.15   9.86   1.01   0.22   5.58 1.00     1884     1828
##    y0[15] 229.64 136.08 396.83 116.95  27.88 668.82 1.00     2036     1881
##    y0[16]   1.48   0.88   1.97   0.74   0.17   4.86 1.00     1685     1574
##    y0[17]   0.15   0.07   0.76   0.07   0.01   0.46 1.00     1649     1821
##    y0[18] 160.10  92.88 253.38  82.21  17.97 478.66 1.00     2110     1989
##    y0[19] 118.89  69.71 224.96  58.31  13.38 351.97 1.00     2084     1868
##    y0[20] 175.03  97.40 310.72  87.05  18.99 555.14 1.00     2126     1611
##    m1[1]   -2.64  -2.64   0.53   0.53  -3.50  -1.75 1.01      872      888
##    m1[2]    1.37   1.37   0.24   0.23   0.98   1.74 1.01      853      998
##    m1[3]    2.54   2.54   0.21   0.20   2.19   2.87 1.00     1200     1212
##    m1[4]    3.24   3.24   0.22   0.20   2.88   3.60 1.00     1654     1396
##    m1[5]    3.75   3.75   0.23   0.22   3.35   4.15 1.00     2051     1487
##    m1[6]    4.15   4.15   0.25   0.24   3.72   4.59 1.00     2268     1624
##    m1[7]    4.48   4.47   0.27   0.25   4.01   4.94 1.00     2327     1596
##    m1[8]    4.75   4.75   0.29   0.27   4.26   5.25 1.00     2309     1651
##    m1[9]    4.99   4.99   0.31   0.29   4.47   5.51 1.00     2275     1596
##    m1[10]   5.21   5.20   0.32   0.31   4.66   5.75 1.00     2233     1546
##    y1[1]    0.13   0.07   0.20   0.06   0.01   0.40 1.00     1434     1436
##    y1[2]    6.36   3.95  10.29   3.36   0.75  18.00 1.00     1753     1853
##    y1[3]   20.34  12.53  27.11  10.50   2.56  60.89 1.00     2158     1935
##    y1[4]   40.41  25.49  53.29  21.06   5.28 121.35 1.00     2118     1894
##    y1[5]   72.14  42.34 104.31  36.22   8.46 227.50 1.00     1965     1810
##    y1[6]  106.71  63.13 143.03  54.91  12.58 327.53 1.00     2081     1965
##    y1[7]  152.74  89.73 356.89  77.32  17.86 448.15 1.00     2249     2055
##    y1[8]  184.17 116.68 237.58  98.08  22.92 562.03 1.00     1731     1863
##    y1[9]  244.00 145.46 311.55 132.50  25.79 781.64 1.00     2225     1930
##    y1[10] 319.77 185.35 611.96 163.76  30.88 922.32 1.00     2012     1930

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

qplot(log(x),log(y),col=I('red'))+
  geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=m),xy,col='black')

exponential increasing/decreasing

y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)  
log y=log b0+b1* x  -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)  
b1>0 x->Infinity,y->Infinity  
b1<0 x->Infinity,y->+0  
n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

ex8-3-1.stan

y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0*exp(b1*x),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*exp(b1*x[i]);
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*exp(b1*x1[i]);
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)
## Initial log joint probability = -1.03578e+08 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##        7      -19.3078   0.000234829   0.000118616           1           1       39    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -19.31
##    b0         0.00
##    b1       -10.00
##    s          1.59
##    m0[1]      0.00
##    m0[2]      0.00
##    m0[3]      0.00
##    m0[4]      0.00
##    m0[5]      0.00
##    m0[6]      0.00
##    m0[7]      0.00
##    m0[8]      0.00
##    m0[9]      0.00
##    m0[10]     0.00
##    m0[11]     0.00
##    m0[12]     0.00
##    m0[13]     0.00
##    m0[14]     0.00
##    m0[15]     0.00
##    m0[16]     0.00
##    m0[17]     0.00
##    m0[18]     0.00
##    m0[19]     0.00
##    m0[20]     0.00
##    y0[1]     -1.06
##    y0[2]      0.29
##    y0[3]     -0.11
##    y0[4]     -0.69
##    y0[5]     -0.80
##    y0[6]     -0.47
##    y0[7]      1.21
##    y0[8]      0.46
##    y0[9]     -1.59
##    y0[10]    -0.52
##    y0[11]    -1.15
##    y0[12]     1.43
##    y0[13]    -0.03
##    y0[14]     0.92
##    y0[15]     1.13
##    y0[16]    -3.54
##    y0[17]    -0.07
##    y0[18]    -0.68
##    y0[19]    -1.11
##    y0[20]    -1.20
##    m1[1]      0.00
##    m1[2]      0.00
##    m1[3]      0.00
##    m1[4]      0.00
##    m1[5]      0.00
##    m1[6]      0.00
##    m1[7]      0.00
##    m1[8]      0.00
##    m1[9]      0.00
##    m1[10]     0.00
##    y1[1]      2.52
##    y1[2]     -0.90
##    y1[3]      1.30
##    y1[4]      2.21
##    y1[5]      1.79
##    y1[6]      1.62
##    y1[7]      0.86
##    y1[8]     -0.35
##    y1[9]      2.94
##    y1[10]     0.58
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__    6.84   7.22 1.40 1.11  4.10  8.35 1.00      456      456
##    b0     14.29  12.96 5.44 3.91  8.32 25.30 1.00      356      318
##    b1     -2.74  -2.66 0.63 0.57 -3.98 -1.85 1.00      366      321
##    s       0.54   0.52 0.10 0.09  0.41  0.71 1.01      491      608
##    m0[1]   0.00   0.00 0.01 0.00  0.00  0.01 1.00      371      355
##    m0[2]   0.00   0.00 0.01 0.00  0.00  0.01 1.00      371      355
##    m0[3]   0.05   0.04 0.05 0.04  0.00  0.15 1.00      377      400
##    m0[4]   0.01   0.01 0.02 0.01  0.00  0.05 1.00      374      379
##    m0[5]   0.00   0.00 0.00 0.00  0.00  0.00 1.00      371      355
##    m0[6]   0.36   0.35 0.17 0.17  0.11  0.68 1.00      392      410
##    m0[7]   0.00   0.00 0.00 0.00  0.00  0.00 1.00      370      333
##    m0[8]   2.08   2.10 0.27 0.25  1.57  2.49 1.00      689      627
##    m0[9]   0.01   0.00 0.01 0.00  0.00  0.03 1.00      373      379
##    m0[10]  0.00   0.00 0.00 0.00  0.00  0.01 1.00      371      355
##    m0[11]  0.02   0.01 0.02 0.01  0.00  0.06 1.00      374      379
##    m0[12]  0.01   0.01 0.02 0.01  0.00  0.05 1.00      374      379
##    m0[13]  4.28   4.25 0.43 0.40  3.60  5.03 1.00      603      794
##    m0[14]  2.31   2.32 0.26 0.24  1.83  2.71 1.00      846      637
##    m0[15]  0.06   0.05 0.06 0.04  0.01  0.18 1.00      377      400
##    m0[16]  1.70   1.73 0.28 0.26  1.18  2.13 1.00      542      552
##    m0[17]  0.00   0.00 0.00 0.00  0.00  0.01 1.00      371      355
##    m0[18]  0.00   0.00 0.00 0.00  0.00  0.00 1.00      370      333
##    m0[19]  3.56   3.56 0.30 0.29  3.10  4.06 1.00     1149     1089
##    m0[20]  1.32   1.33 0.28 0.27  0.81  1.77 1.00      468      521
##    y0[1]  -0.02  -0.01 0.55 0.53 -0.92  0.88 1.00     2065     2010
##    y0[2]   0.01   0.02 0.54 0.51 -0.88  0.87 1.00     2090     1782
##    y0[3]   0.05   0.05 0.55 0.53 -0.86  0.93 1.00     2033     1848
##    y0[4]   0.01  -0.02 0.55 0.53 -0.86  0.91 1.00     1959     1857
##    y0[5]  -0.01   0.00 0.56 0.53 -0.94  0.88 1.00     2117     1972
##    y0[6]   0.37   0.35 0.59 0.57 -0.61  1.31 1.00     1583     1738
##    y0[7]   0.03   0.03 0.55 0.52 -0.88  0.92 1.00     1763     1703
##    y0[8]   2.07   2.08 0.62 0.58  1.01  3.04 1.00     1625     1599
##    y0[9]   0.03   0.01 0.55 0.51 -0.86  0.92 1.00     1958     1678
##    y0[10] -0.01  -0.01 0.54 0.52 -0.90  0.87 1.00     2007     1883
##    y0[11]  0.01   0.02 0.55 0.52 -0.85  0.91 1.00     2014     1892
##    y0[12]  0.01   0.03 0.54 0.53 -0.87  0.85 1.00     1804     1694
##    y0[13]  4.28   4.27 0.70 0.69  3.18  5.49 1.00     1058     1088
##    y0[14]  2.30   2.33 0.61 0.56  1.25  3.30 1.00     1289     1564
##    y0[15]  0.06   0.06 0.52 0.51 -0.78  0.90 1.00     1910     1892
##    y0[16]  1.70   1.72 0.61 0.56  0.70  2.67 1.00     1296     1528
##    y0[17]  0.01   0.01 0.54 0.51 -0.88  0.88 1.00     1986     1877
##    y0[18]  0.01   0.01 0.54 0.51 -0.87  0.92 1.00     1976     1973
##    y0[19]  3.55   3.56 0.61 0.58  2.53  4.55 1.00     1454     1306
##    y0[20]  1.33   1.34 0.62 0.59  0.32  2.30 1.00     1420     1595
##    m1[1]   4.28   4.25 0.43 0.40  3.60  5.03 1.00      603      794
##    m1[2]   1.15   1.16 0.27 0.27  0.66  1.59 1.00      447      454
##    m1[3]   0.33   0.32 0.16 0.16  0.10  0.64 1.00      390      410
##    m1[4]   0.10   0.09 0.08 0.06  0.01  0.26 1.00      380      400
##    m1[5]   0.03   0.02 0.04 0.02  0.00  0.11 1.00      375      379
##    m1[6]   0.01   0.01 0.02 0.01  0.00  0.04 1.00      373      379
##    m1[7]   0.00   0.00 0.01 0.00  0.00  0.02 1.00      372      355
##    m1[8]   0.00   0.00 0.00 0.00  0.00  0.01 1.00      371      355
##    m1[9]   0.00   0.00 0.00 0.00  0.00  0.00 1.00      370      333
##    m1[10]  0.00   0.00 0.00 0.00  0.00  0.00 1.00      370      333
##    y1[1]   4.27   4.26 0.68 0.66  3.17  5.40 1.00     1211     1220
##    y1[2]   1.18   1.18 0.62 0.59  0.13  2.15 1.00     1180     1120
##    y1[3]   0.33   0.32 0.57 0.55 -0.58  1.28 1.00     1563     1505
##    y1[4]   0.09   0.10 0.56 0.54 -0.83  1.00 1.00     2113     1447
##    y1[5]   0.04   0.03 0.54 0.52 -0.84  0.93 1.00     2012     1794
##    y1[6]   0.03   0.04 0.53 0.53 -0.81  0.88 1.00     2085     1892
##    y1[7]   0.00   0.00 0.54 0.51 -0.91  0.88 1.00     2073     1706
##    y1[8]  -0.02  -0.01 0.55 0.55 -0.92  0.87 1.00     2030     1882
##    y1[9]   0.00   0.00 0.54 0.51 -0.86  0.87 1.00     2009     1887
##    y1[10] -0.01  -0.02 0.54 0.52 -0.89  0.88 1.00     2084     1684

ex8-3-2.stan

log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real<lower=-10,upper=10> b1;
  real<lower=0> s;
}
model{
  y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=log(b0)+b1*x[i];
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=log(b0)+b1*x1[i];
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -1167.44 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       36       -14.706   0.000243316   0.000413939           1           1       56    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -14.71
##    b0        10.50
##    b1        -2.03
##    s          0.50
##    m0[1]     -6.99
##    m0[2]      0.60
##    m0[3]     -4.02
##    m0[4]     -6.57
##    m0[5]     -3.17
##    m0[6]      0.44
##    m0[7]     -4.85
##    m0[8]     -0.78
##    m0[9]      1.80
##    m0[10]    -1.18
##    m0[11]    -5.96
##    m0[12]     1.56
##    m0[13]    -6.64
##    m0[14]    -0.34
##    m0[15]    -6.65
##    m0[16]    -3.46
##    m0[17]    -1.27
##    m0[18]     1.37
##    m0[19]    -5.78
##    m0[20]    -3.70
##    y0[1]      0.00
##    y0[2]      2.44
##    y0[3]      0.03
##    y0[4]      0.00
##    y0[5]      0.05
##    y0[6]      1.09
##    y0[7]      0.01
##    y0[8]      0.54
##    y0[9]      3.55
##    y0[10]     0.43
##    y0[11]     0.00
##    y0[12]     7.80
##    y0[13]     0.00
##    y0[14]     0.95
##    y0[15]     0.00
##    y0[16]     0.04
##    y0[17]     0.24
##    y0[18]     3.79
##    y0[19]     0.00
##    y0[20]     0.01
##    m1[1]      1.80
##    m1[2]      0.82
##    m1[3]     -0.16
##    m1[4]     -1.13
##    m1[5]     -2.11
##    m1[6]     -3.08
##    m1[7]     -4.06
##    m1[8]     -5.04
##    m1[9]     -6.01
##    m1[10]    -6.99
##    y1[1]      9.00
##    y1[2]      1.75
##    y1[3]      0.73
##    y1[4]      0.21
##    y1[5]      0.13
##    y1[6]      0.02
##    y1[7]      0.04
##    y1[8]      0.00
##    y1[9]      0.00
##    y1[10]     0.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -13.12 -12.80 1.34 1.16 -15.69 -11.63 1.00      632      894
##    b0      11.42  10.94 3.11 2.68   7.38  17.27 1.01      483      529
##    b1      -2.04  -2.04 0.09 0.08  -2.19  -1.90 1.00      639      818
##    s        0.58   0.57 0.11 0.10   0.43   0.77 1.01      650      840
##    m0[1]   -7.01  -7.00 0.23 0.22  -7.39  -6.64 1.00     1657     1535
##    m0[2]    0.64   0.62 0.20 0.19   0.32   0.99 1.01      502      554
##    m0[3]   -4.02  -4.02 0.15 0.14  -4.27  -3.78 1.00     2002     1261
##    m0[4]   -6.58  -6.58 0.21 0.20  -6.94  -6.24 1.00     1776     1614
##    m0[5]   -3.17  -3.17 0.14 0.13  -3.39  -2.93 1.00     1380     1189
##    m0[6]    0.48   0.47 0.19 0.18   0.17   0.83 1.01      505      570
##    m0[7]   -4.85  -4.85 0.17 0.15  -5.13  -4.59 1.00     2212     1517
##    m0[8]   -0.76  -0.76 0.16 0.15  -1.00  -0.48 1.00      571      623
##    m0[9]    1.84   1.83 0.24 0.23   1.47   2.26 1.01      486      518
##    m0[10]  -1.16  -1.16 0.15 0.14  -1.40  -0.89 1.00      616      709
##    m0[11]  -5.97  -5.97 0.20 0.18  -6.30  -5.65 1.00     1968     1616
##    m0[12]   1.60   1.59 0.23 0.22   1.24   2.00 1.01      487      518
##    m0[13]  -6.66  -6.65 0.22 0.20  -7.02  -6.30 1.00     1755     1614
##    m0[14]  -0.30  -0.31 0.17 0.16  -0.57   0.00 1.00      538      636
##    m0[15]  -6.67  -6.66 0.22 0.20  -7.02  -6.31 1.00     1753     1614
##    m0[16]  -3.46  -3.46 0.14 0.13  -3.69  -3.23 1.00     1601     1269
##    m0[17]  -1.25  -1.26 0.15 0.14  -1.49  -0.98 1.00      629      709
##    m0[18]   1.41   1.40 0.22 0.21   1.05   1.81 1.01      488      518
##    m0[19]  -5.79  -5.79 0.19 0.18  -6.11  -5.49 1.00     2024     1589
##    m0[20]  -3.69  -3.69 0.14 0.13  -3.93  -3.46 1.00     1777     1284
##    y0[1]    0.00   0.00 0.00 0.00   0.00   0.00 1.00     2064     1895
##    y0[2]    2.29   1.84 1.75 1.02   0.69   5.27 1.00     1316     1799
##    y0[3]    0.02   0.02 0.01 0.01   0.01   0.05 1.00     1987     1843
##    y0[4]    0.00   0.00 0.00 0.00   0.00   0.00 1.00     1957     1792
##    y0[5]    0.05   0.04 0.03 0.02   0.02   0.11 1.00     2017     1879
##    y0[6]    2.01   1.67 1.49 0.94   0.62   4.52 1.00     1611     1799
##    y0[7]    0.01   0.01 0.01 0.00   0.00   0.02 1.00     1931     1919
##    y0[8]    0.58   0.48 0.41 0.26   0.18   1.35 1.00     1646     1893
##    y0[9]    7.85   6.32 5.81 3.65   2.31  18.10 1.00     1378     1095
##    y0[10]   0.38   0.31 0.27 0.18   0.11   0.86 1.00     1789     1680
##    y0[11]   0.00   0.00 0.00 0.00   0.00   0.01 1.00     1806     1933
##    y0[12]   5.97   4.88 4.24 2.69   1.79  13.91 1.00     1445     1516
##    y0[13]   0.00   0.00 0.00 0.00   0.00   0.00 1.00     1986     1514
##    y0[14]   0.91   0.74 0.79 0.43   0.26   2.05 1.00     1748     1448
##    y0[15]   0.00   0.00 0.00 0.00   0.00   0.00 1.00     2054     1773
##    y0[16]   0.04   0.03 0.03 0.02   0.01   0.08 1.00     2130     1838
##    y0[17]   0.34   0.28 0.24 0.16   0.10   0.76 1.00     1624     1919
##    y0[18]   5.07   4.12 3.86 2.35   1.52  11.78 1.00     1360     1526
##    y0[19]   0.00   0.00 0.00 0.00   0.00   0.01 1.00     2051     1816
##    y0[20]   0.03   0.03 0.02 0.01   0.01   0.07 1.00     1748     1915
##    m1[1]    1.84   1.83 0.24 0.23   1.47   2.26 1.01      486      518
##    m1[2]    0.86   0.85 0.21 0.19   0.53   1.22 1.01      497      541
##    m1[3]   -0.12  -0.13 0.18 0.17  -0.40   0.19 1.00      528      636
##    m1[4]   -1.11  -1.11 0.16 0.14  -1.35  -0.84 1.00      610      683
##    m1[5]   -2.09  -2.09 0.14 0.13  -2.31  -1.85 1.00      816      906
##    m1[6]   -3.07  -3.08 0.14 0.13  -3.30  -2.84 1.00     1313     1179
##    m1[7]   -4.06  -4.06 0.15 0.14  -4.31  -3.82 1.00     1999     1237
##    m1[8]   -5.04  -5.04 0.17 0.16  -5.33  -4.77 1.00     2203     1513
##    m1[9]   -6.02  -6.02 0.20 0.18  -6.35  -5.71 1.00     1953     1616
##    m1[10]  -7.01  -7.00 0.23 0.22  -7.39  -6.64 1.00     1657     1535
##    y1[1]    7.80   6.32 5.96 3.46   2.27  18.72 1.00     1376     1585
##    y1[2]    2.80   2.32 1.92 1.33   0.83   6.54 1.00     1575     1681
##    y1[3]    1.06   0.87 0.86 0.49   0.31   2.32 1.00     1654     1506
##    y1[4]    0.40   0.33 0.26 0.18   0.13   0.88 1.00     1846     1951
##    y1[5]    0.15   0.12 0.10 0.07   0.05   0.34 1.00     1801     1891
##    y1[6]    0.05   0.05 0.04 0.03   0.02   0.12 1.00     2081     1841
##    y1[7]    0.02   0.02 0.02 0.01   0.01   0.05 1.00     2264     1866
##    y1[8]    0.01   0.01 0.01 0.00   0.00   0.02 1.00     1725     1929
##    y1[9]    0.00   0.00 0.00 0.00   0.00   0.01 1.00     1950     1732
##    y1[10]   0.00   0.00 0.00 0.00   0.00   0.00 1.00     1933     1691

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,log(y),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

growth curve

y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)  
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)

ex8-4.stan

growth curve

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=10> b1;
  real<lower=0> s;
}
model{
  y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(1-exp(-b1*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(1-exp(-b1*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)
## Initial log joint probability = -3396.36 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -49.0744    0.00344815    0.00157874           1           1       15    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -49.07
##    b0         0.00
##    b1         1.84
##    s          7.05
##    m0[1]      0.00
##    m0[2]      0.00
##    m0[3]      0.00
##    m0[4]      0.00
##    m0[5]      0.00
##    m0[6]      0.00
##    m0[7]      0.00
##    m0[8]      0.00
##    m0[9]      0.00
##    m0[10]     0.00
##    m0[11]     0.00
##    m0[12]     0.00
##    m0[13]     0.00
##    m0[14]     0.00
##    m0[15]     0.00
##    m0[16]     0.00
##    m0[17]     0.00
##    m0[18]     0.00
##    m0[19]     0.00
##    m0[20]     0.00
##    y0[1]     -8.04
##    y0[2]     10.33
##    y0[3]      3.92
##    y0[4]      4.07
##    y0[5]     -4.23
##    y0[6]      4.83
##    y0[7]     11.15
##    y0[8]     -1.18
##    y0[9]     -3.67
##    y0[10]    -0.31
##    y0[11]     5.87
##    y0[12]   -17.71
##    y0[13]    -4.30
##    y0[14]    -0.12
##    y0[15]   -12.58
##    y0[16]     8.36
##    y0[17]     1.83
##    y0[18]     6.85
##    y0[19]     3.86
##    y0[20]     2.79
##    m1[1]      0.00
##    m1[2]      0.00
##    m1[3]      0.00
##    m1[4]      0.00
##    m1[5]      0.00
##    m1[6]      0.00
##    m1[7]      0.00
##    m1[8]      0.00
##    m1[9]      0.00
##    m1[10]     0.00
##    y1[1]     -4.59
##    y1[2]    -10.28
##    y1[3]    -10.30
##    y1[4]      0.40
##    y1[5]      0.81
##    y1[6]     -1.43
##    y1[7]     -4.01
##    y1[8]     -3.04
##    y1[9]      1.49
##    y1[10]     3.87
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -11.47 -11.10 1.37 1.05 -14.17 -10.03 1.00      642      798
##    b0       9.70   9.60 0.86 0.75   8.51  11.26 1.00      737      783
##    b1       0.44   0.43 0.10 0.09   0.29   0.61 1.01      621      682
##    s        1.20   1.17 0.22 0.19   0.90   1.61 1.00      986      839
##    m0[1]    1.60   1.58 0.22 0.21   1.25   1.99 1.00      626      760
##    m0[2]    4.78   4.76 0.44 0.42   4.06   5.52 1.00      666      853
##    m0[3]    6.54   6.54 0.40 0.39   5.88   7.18 1.00      838      898
##    m0[4]    5.05   5.03 0.44 0.42   4.33   5.78 1.00      675      881
##    m0[5]    6.69   6.69 0.39 0.38   6.04   7.31 1.00      882      915
##    m0[6]    6.74   6.74 0.39 0.38   6.09   7.35 1.00      897      971
##    m0[7]    9.16   9.17 0.51 0.49   8.33  10.01 1.00      961     1204
##    m0[8]    7.26   7.27 0.35 0.34   6.68   7.83 1.00     1211     1015
##    m0[9]    9.05   9.06 0.47 0.45   8.28   9.82 1.00     1052     1286
##    m0[10]   9.16   9.17 0.51 0.49   8.33  10.02 1.00      958     1204
##    m0[11]   4.72   4.70 0.44 0.42   4.00   5.45 1.00      665      852
##    m0[12]   8.09   8.08 0.32 0.31   7.58   8.62 1.00     1990     1195
##    m0[13]   8.87   8.88 0.42 0.39   8.19   9.55 1.00     1239     1161
##    m0[14]   4.74   4.72 0.44 0.42   4.02   5.48 1.00      665      852
##    m0[15]   7.69   7.69 0.33 0.31   7.16   8.23 1.00     1659     1174
##    m0[16]   9.19   9.19 0.52 0.51   8.35  10.06 1.00      941     1204
##    m0[17]   1.83   1.81 0.25 0.24   1.44   2.27 1.00      627      779
##    m0[18]   8.97   8.98 0.45 0.43   8.25   9.71 1.00     1126     1269
##    m0[19]   8.25   8.24 0.33 0.31   7.73   8.78 1.00     1999     1184
##    m0[20]   1.19   1.17 0.17 0.16   0.92   1.49 1.00      625      759
##    y0[1]    1.59   1.57 1.21 1.17  -0.36   3.60 1.00     1956     1645
##    y0[2]    4.73   4.74 1.31 1.26   2.57   6.88 1.00     1757     1997
##    y0[3]    6.52   6.54 1.30 1.24   4.39   8.70 1.00     2007     2012
##    y0[4]    5.05   5.04 1.30 1.26   2.96   7.26 1.00     1862     1773
##    y0[5]    6.67   6.64 1.26 1.21   4.68   8.74 1.00     1993     1880
##    y0[6]    6.73   6.72 1.30 1.24   4.57   8.92 1.00     1581     1892
##    y0[7]    9.15   9.14 1.28 1.26   7.06  11.29 1.00     1751     1777
##    y0[8]    7.25   7.25 1.26 1.21   5.15   9.34 1.00     1800     1932
##    y0[9]    9.07   9.10 1.31 1.28   6.89  11.19 1.00     1994     1903
##    y0[10]   9.16   9.12 1.32 1.24   7.01  11.41 1.00     1647     1808
##    y0[11]   4.74   4.74 1.27 1.22   2.68   6.89 1.00     1522     1711
##    y0[12]   8.08   8.06 1.29 1.16   6.03  10.23 1.00     1651     1821
##    y0[13]   8.84   8.85 1.31 1.24   6.67  10.98 1.00     1881     1657
##    y0[14]   4.75   4.75 1.31 1.24   2.60   6.87 1.00     1434     1644
##    y0[15]   7.67   7.70 1.30 1.26   5.60   9.74 1.00     1856     1734
##    y0[16]   9.17   9.18 1.31 1.28   6.98  11.32 1.00     1865     1822
##    y0[17]   1.84   1.81 1.25 1.20  -0.18   3.86 1.00     1971     1820
##    y0[18]   8.96   9.00 1.28 1.23   6.77  11.01 1.00     1691     1840
##    y0[19]   8.28   8.30 1.28 1.21   6.15  10.38 1.00     2064     2099
##    y0[20]   1.20   1.17 1.25 1.20  -0.78   3.18 1.00     2112     1969
##    m1[1]    1.19   1.17 0.17 0.16   0.92   1.49 1.00      625      759
##    m1[2]    3.60   3.58 0.40 0.38   2.97   4.30 1.00      642      793
##    m1[3]    5.31   5.29 0.44 0.43   4.58   6.04 1.00      687      881
##    m1[4]    6.52   6.52 0.40 0.39   5.85   7.16 1.00      835      898
##    m1[5]    7.38   7.39 0.35 0.33   6.83   7.94 1.00     1330     1134
##    m1[6]    8.00   7.99 0.32 0.30   7.48   8.52 1.00     1942     1222
##    m1[7]    8.45   8.45 0.35 0.32   7.90   9.00 1.00     1840     1253
##    m1[8]    8.78   8.79 0.40 0.38   8.13   9.42 1.00     1354     1287
##    m1[9]    9.01   9.02 0.46 0.44   8.27   9.76 1.00     1084     1268
##    m1[10]   9.19   9.19 0.52 0.51   8.35  10.06 1.00      941     1204
##    y1[1]    1.21   1.21 1.25 1.17  -0.83   3.26 1.00     1711     1840
##    y1[2]    3.59   3.60 1.29 1.25   1.42   5.68 1.00     1724     1893
##    y1[3]    5.32   5.33 1.31 1.25   3.21   7.45 1.00     1666     1590
##    y1[4]    6.52   6.49 1.33 1.25   4.38   8.72 1.00     1618     1747
##    y1[5]    7.44   7.44 1.26 1.29   5.36   9.45 1.00     1707     1998
##    y1[6]    7.98   7.99 1.28 1.22   5.81  10.02 1.00     1959     1906
##    y1[7]    8.48   8.48 1.30 1.29   6.30  10.55 1.00     1646     1353
##    y1[8]    8.80   8.81 1.29 1.25   6.70  10.92 1.00     2038     1972
##    y1[9]    9.04   9.03 1.28 1.23   6.91  11.11 1.00     1945     1959
##    y1[10]   9.19   9.19 1.32 1.27   7.09  11.32 1.00     1607     1314

sigmoid curve

y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)

ex8-5.stan

sigmoid curve

data{
  int N;
  vector[N] x;
  int Ym;
  array[N] int y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=100> b1;
}
model{
  y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
  array[N] int y0;
  for(i in 1:N){
    y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
  }
  array[N1] int y1;
  for(i in 1:N1){
    y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
  }
}
n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-5.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "y0"   "y1"
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -59.64 -59.26 1.13 0.75 -62.00 -58.61 1.00      641      796
##    b0       3.92   3.93 0.15 0.15   3.66   4.17 1.00     2268     1417
##    b1       1.58   1.57 0.24 0.23   1.21   1.99 1.00      564      491
##    y0[1]    9.80  10.00 0.46 0.00   9.00  10.00 1.00     1847       NA
##    y0[2]    9.98  10.00 0.14 0.00  10.00  10.00 1.00     1781       NA
##    y0[3]    9.57  10.00 0.68 0.00   8.00  10.00 1.00     1438       NA
##    y0[4]    0.06   0.00 0.26 0.00   0.00   1.00 1.00     1645     1656
##    y0[5]    7.86   8.00 1.34 1.48   6.00  10.00 1.00     1179       NA
##    y0[6]    0.06   0.00 0.24 0.00   0.00   1.00 1.00     1521     1504
##    y0[7]    9.83  10.00 0.43 0.00   9.00  10.00 1.00     1868       NA
##    y0[8]    9.80  10.00 0.45 0.00   9.00  10.00 1.00     1765       NA
##    y0[9]    0.12   0.00 0.36 0.00   0.00   1.00 1.00     1659     1646
##    y0[10]   2.47   2.00 1.47 1.48   0.00   5.00 1.00     1882     1693
##    y0[11]   1.15   1.00 1.05 1.48   0.00   3.00 1.00     1790     1978
##    y0[12]   2.37   2.00 1.43 1.48   0.00   5.00 1.00     1795     1914
##    y0[13]   0.21   0.00 0.47 0.00   0.00   1.00 1.00     1931     1838
##    y0[14]   4.97   5.00 1.65 1.48   2.00   8.00 1.00     1958     1776
##    y0[15]   8.50   9.00 1.20 1.48   6.00  10.00 1.00     1423       NA
##    y0[16]   4.58   5.00 1.62 1.48   2.00   7.00 1.00     2077     1813
##    y0[17]   0.69   0.00 0.83 0.00   0.00   2.00 1.00     1862     1705
##    y0[18]   5.73   6.00 1.64 1.48   3.00   8.00 1.00     1948     2019
##    y0[19]   1.91   2.00 1.30 1.48   0.00   4.00 1.00     1860     1927
##    y0[20]   9.37  10.00 0.81 0.00   8.00  10.00 1.00     1863       NA
##    y1[1]    0.06   0.00 0.25 0.00   0.00   1.00 1.00     1666     1440
##    y1[2]    0.17   0.00 0.42 0.00   0.00   1.00 1.00     1964     1998
##    y1[3]    0.63   0.00 0.79 0.00   0.00   2.00 1.00     1454     1506
##    y1[4]    2.05   2.00 1.37 1.48   0.00   5.00 1.00     1802     1841
##    y1[5]    4.91   5.00 1.68 1.48   2.00   8.00 1.00     2051     2085
##    y1[6]    7.89   8.00 1.37 1.48   5.00  10.00 1.00     1831       NA
##    y1[7]    9.33  10.00 0.87 0.00   8.00  10.00 1.00     1490       NA
##    y1[8]    9.80  10.00 0.45 0.00   9.00  10.00 1.00     1886       NA
##    y1[9]    9.95  10.00 0.23 0.00   9.00  10.00 1.00     1963       NA
##    y1[10]   9.99  10.00 0.12 0.00  10.00  10.00 1.00     1935       NA

y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=ymed),xy,col='black')

up and down

y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))

ex8-6.stan

up and down

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=200> b0;
  real<lower=0,upper=1> b1;
  real<lower=0,upper=1> b2;
  real<lower=0,upper=10> s;
}
model{
  y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-6.stan')

fn(mdl,data)
## Initial log joint probability = -574.473 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1eEBAt/model-7857f9e47b8.stan', line 15, column 2 to column 41) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       55      -5.86343   2.19592e-05     0.0029161           1           1      137    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -5.86
##    b0        98.77
##    b1         0.03
##    b2         0.20
##    s          0.81
##    m0[1]     35.51
##    m0[2]     48.72
##    m0[3]     34.66
##    m0[4]     25.66
##    m0[5]     14.61
##    m0[6]     37.52
##    m0[7]     58.04
##    m0[8]     58.46
##    m0[9]     25.17
##    m0[10]    23.49
##    m0[11]    58.10
##    m0[12]    60.53
##    m0[13]    55.61
##    m0[14]    26.04
##    m0[15]    53.35
##    m0[16]    38.28
##    m0[17]    57.80
##    m0[18]    23.54
##    m0[19]    56.38
##    m0[20]    46.36
##    y0[1]     36.37
##    y0[2]     48.25
##    y0[3]     35.77
##    y0[4]     26.26
##    y0[5]     15.61
##    y0[6]     38.06
##    y0[7]     59.13
##    y0[8]     57.54
##    y0[9]     25.62
##    y0[10]    22.49
##    y0[11]    57.84
##    y0[12]    59.27
##    y0[13]    53.97
##    y0[14]    26.20
##    y0[15]    53.73
##    y0[16]    38.44
##    y0[17]    58.95
##    y0[18]    24.36
##    y0[19]    57.09
##    y0[20]    45.07
##    m1[1]     14.61
##    m1[2]     54.50
##    m1[3]     60.86
##    m1[4]     56.89
##    m1[5]     50.29
##    m1[6]     43.57
##    m1[7]     37.45
##    m1[8]     32.09
##    m1[9]     27.46
##    m1[10]    23.49
##    y1[1]     14.56
##    y1[2]     54.25
##    y1[3]     61.37
##    y1[4]     56.19
##    y1[5]     50.03
##    y1[6]     44.03
##    y1[7]     36.82
##    y1[8]     31.69
##    y1[9]     26.92
##    y1[10]    24.38
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable  mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -9.81  -9.44 1.58 1.36 -12.72  -7.97 1.00      443      565
##    b0     98.77  98.73 2.35 2.18  94.93 102.74 1.00      919      550
##    b1      0.03   0.03 0.00 0.00   0.03   0.03 1.00      943      692
##    b2      0.20   0.20 0.01 0.01   0.19   0.22 1.00      966      662
##    s       0.95   0.92 0.18 0.16   0.72   1.31 1.02      256      304
##    m0[1]  35.52  35.53 0.28 0.27  35.06  35.99 1.00     1847     1602
##    m0[2]  48.74  48.73 0.63 0.61  47.76  49.80 1.00     1245     1317
##    m0[3]  34.67  34.67 0.29 0.27  34.20  35.14 1.00     1816     1454
##    m0[4]  25.68  25.68 0.35 0.33  25.09  26.26 1.00     1269     1173
##    m0[5]  14.62  14.61 0.33 0.32  14.11  15.19 1.00     1060      970
##    m0[6]  37.53  37.54 0.28 0.27  37.07  37.99 1.00     1899     1620
##    m0[7]  58.04  58.06 0.36 0.34  57.44  58.59 1.00     1264     1081
##    m0[8]  58.46  58.48 0.36 0.34  57.87  59.00 1.00     1282     1057
##    m0[9]  25.18  25.18 0.35 0.33  24.59  25.77 1.00     1252     1159
##    m0[10] 23.51  23.51 0.36 0.34  22.90  24.12 1.00     1210     1028
##    m0[11] 58.10  58.10 0.48 0.47  57.35  58.92 1.00     1875     1354
##    m0[12] 60.52  60.54 0.35 0.33  59.93  61.06 1.00     1550     1366
##    m0[13] 55.60  55.62 0.36 0.34  55.03  56.16 1.00     1222     1080
##    m0[14] 26.06  26.06 0.35 0.33  25.48  26.64 1.00     1282     1174
##    m0[15] 53.35  53.37 0.35 0.33  52.78  53.90 1.00     1243     1082
##    m0[16] 38.29  38.29 0.28 0.27  37.83  38.76 1.00     1906     1593
##    m0[17] 57.81  57.81 0.49 0.48  57.04  58.64 1.00     1861     1393
##    m0[18] 23.56  23.56 0.36 0.34  22.95  24.17 1.00     1212     1044
##    m0[19] 56.37  56.39 0.36 0.34  55.80  56.92 1.00     1228     1092
##    m0[20] 46.36  46.37 0.31 0.28  45.84  46.85 1.00     1511     1282
##    y0[1]  35.55  35.56 1.02 0.97  33.96  37.24 1.00     1908     1858
##    y0[2]  48.72  48.76 1.17 1.11  46.81  50.63 1.00     1801     1776
##    y0[3]  34.67  34.70 1.01 0.96  32.97  36.31 1.00     1686     1642
##    y0[4]  25.67  25.67 1.02 0.97  24.01  27.35 1.00     1636     1582
##    y0[5]  14.60  14.60 1.03 1.02  12.96  16.28 1.00     1876     1704
##    y0[6]  37.55  37.55 1.01 0.96  35.90  39.22 1.00     1905     1919
##    y0[7]  58.04  58.00 1.05 1.03  56.36  59.72 1.00     1830     1518
##    y0[8]  58.48  58.51 1.05 1.00  56.76  60.21 1.00     1853     1555
##    y0[9]  25.19  25.21 1.02 0.95  23.49  26.88 1.00     1948     1733
##    y0[10] 23.51  23.50 1.03 0.97  21.83  25.16 1.00     1931     1912
##    y0[11] 58.12  58.12 1.10 1.02  56.34  59.94 1.00     1703     1808
##    y0[12] 60.52  60.53 1.02 0.96  58.82  62.17 1.00     1951     1572
##    y0[13] 55.61  55.61 1.06 0.99  53.92  57.37 1.00     1610     1869
##    y0[14] 26.04  26.06 1.04 0.98  24.26  27.76 1.00     1974     1729
##    y0[15] 53.37  53.37 1.05 0.98  51.70  55.10 1.00     1895     1724
##    y0[16] 38.28  38.31 1.04 0.99  36.56  39.90 1.00     2085     1932
##    y0[17] 57.78  57.76 1.09 1.06  55.96  59.50 1.00     1950     1739
##    y0[18] 23.55  23.53 1.03 1.00  21.88  25.25 1.00     1998     1874
##    y0[19] 56.34  56.31 1.04 0.99  54.68  58.07 1.00     1820     1907
##    y0[20] 46.35  46.34 1.01 0.97  44.71  48.00 1.00     2046     2057
##    m1[1]  14.62  14.61 0.33 0.32  14.11  15.19 1.00     1060      970
##    m1[2]  54.51  54.51 0.56 0.54  53.61  55.44 1.00     1542     1378
##    m1[3]  60.86  60.87 0.36 0.34  60.26  61.41 1.00     1729     1474
##    m1[4]  56.88  56.90 0.36 0.34  56.29  57.43 1.00     1235     1092
##    m1[5]  50.29  50.30 0.33 0.32  49.76  50.80 1.00     1319     1271
##    m1[6]  43.57  43.57 0.29 0.27  43.08  44.04 1.00     1721     1362
##    m1[7]  37.46  37.46 0.28 0.27  36.99  37.92 1.00     1898     1707
##    m1[8]  32.10  32.10 0.30 0.28  31.61  32.60 1.00     1603     1354
##    m1[9]  27.48  27.48 0.33 0.32  26.92  28.04 1.00     1335     1135
##    m1[10] 23.51  23.51 0.36 0.34  22.90  24.12 1.00     1210     1028
##    y1[1]  14.61  14.60 1.03 1.00  12.93  16.31 1.00     2009     1645
##    y1[2]  54.47  54.50 1.13 1.05  52.57  56.21 1.00     1893     1370
##    y1[3]  60.86  60.85 1.06 1.01  59.17  62.61 1.00     2069     1892
##    y1[4]  56.89  56.93 1.02 0.99  55.24  58.54 1.00     2006     1823
##    y1[5]  50.27  50.26 1.05 1.00  48.52  51.98 1.00     1817     1702
##    y1[6]  43.56  43.57 1.01 0.96  41.90  45.20 1.00     1801     1653
##    y1[7]  37.47  37.47 1.00 0.94  35.83  39.12 1.00     2123     1770
##    y1[8]  32.12  32.12 0.97 0.94  30.55  33.68 1.00     2093     1810
##    y1[9]  27.47  27.48 1.02 0.95  25.81  29.16 1.00     1922     1740
##    y1[10] 23.48  23.46 1.06 1.04  21.76  25.17 1.00     1491     1530